%dXt = sin(t)dW 
clear all
% index J for X
alphaK=15; alphas=[zeros(1,alphaK);eye(alphaK)];
% time T
T = 1/2;
% time mesh 
num=50; tspan = linspace(0,T,num+1);
%sigma
s=mihandle;
s.sigma=@(s) sin(s);
%options = odeset('AbsTol',1e-12,'RelTol',1e-6);
options = [];
row=size(alphas,1); X0=zeros(row,1);
% solver of X
disp('solving X')

[time,X] = ode45(@randomsg,tspan,X0,options,T,alphaK,s);

hold on
Npath=100000;
r=0.1; K=1;

%%% old code, test variance
% Y=X.*X;  tmp=zeros(size(Y(:,1)));
% for ii=1:size(Y,2)
%     tmp=tmp+Y(:,ii);
% end
% 
% plot(time,tmp,'r-');

%%% E|X|
for jj=1:51
%Y=X(31,:);
Y=X(jj,:);
Z=Y(1)*ones(Npath,1);
for ii=2:16
Z=Z+Y(ii)*randn(Npath,1);
end
price(jj)=sum(abs(Z))/Npath;
end

% discount = exp(-r*time);
% for jj=1:51
% %Y=X(31,:);
% Y=X(jj,:);
% Z=Y(1)*ones(Npath,1);
% for ii=2:16
% Z=Z+Y(ii)*randn(Npath,1);
% end
% Z=exp(Z)-K;
% price(jj)=discount(jj)*sum(Z(Z>0))/Npath;
% end


plot(time,price)


myfun=@(s) sin(s);
paths=zeros(Npath,51);
dt=T/5000;
for ii=2:51
    curtime=0.01*(ii-2)+(1:100)*dt;
    cursigma=myfun(curtime);
    paths(:,ii)=paths(:,ii-1)+(randn(Npath,100)*cursigma')*sqrt(dt);    
end

%%% MC for E|X|
for ii=1:51
    %mcvar=sum(abs(paths(:,31)),1)/Npath
    mcvar(ii)=sum(abs(paths(:,ii)),1)/Npath;
end

% for ii=1:51
%     %mcvar=sum(abs(paths(:,31)),1)/Npath
%     Z=exp(paths(:,ii))-K;
%     mcvar(ii)=discount(ii)*sum(Z(Z>0),1)/Npath;
% end
    plot((0:50)*0.5/50,mcvar,'r')
